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1 Introduction 

There is an important class of computer programs that do calculations in quan- 
tum chromodynamics (QCD) in which the calculation is performed at next-to-leading 
order in perturbation theory and allows for the determination of a variety of char- 
acteristics of the final state. This talk reviews a program of this class in which a 
"completely numerical" integration algorithm is used. 

I consider the calculation of "three-jet-like" observables in e^e~ annihilation. A 
program that does this can be used to calculate a jet cross section (with any in- 
frared safe choice of jet definition) or observables like the thrust distribution. Such 
a program generates random partonic events consisting of three or four final state 
quarks, antiquarks, and gluons. Each event comes with a calculated weight. A sepa- 
rate routine then calculates the contribution to the desired observable for each event, 
averaging over the events with their weights. 

The weights are treated as probabilities. However, these weights can be both 
positive or negative. This is an almost inevitable consequence of quantum mechanics. 
The calculated observable is proportional to the square of a quantum amplitude 
and is thus positive. However, as soon as one divides the amplitude into pieces for 
purposes of calculation, one finds that, while the square of each piece is positive, 
the interference terms between different pieces can have either sign. Thus the kind 
of program discussed here stands in contrast to the tree-level event generators in 
which, by simplifying the physics, one can generally arrange to have all the weights 
be positive, or, even, all be equal to 1. 

To understand the algorithms used in the class of programs described above, it is 
best to think of the calculations as performing integrations over momenta in which the 
quantum matrix elements and the measurement functions form the integrand. There 
are two basic algorithms for performing the integrations. The older is due to Ellis, 
Ross, and Terrano (ERT) |T| . In this method, some of the integrations are performed 
analytically ahead of time. The other integrations are performed numerically by the 
Monte Carlo method. The integrations are divergent and are regulated by analytical 
continuation to 3 — 2e space dimensions and a scheme of subtractions or cutoffs. 
The second method is much newer 0,^,0. In this method, all of the momentum 
integrations are done by Monte Carlo numerical integration. With this method, 
the integrals are all convergent (after removal of the ultraviolet divergences by a 
straightforward renormalization procedure). 

In its current incarnation, the numerical method is not as good as older programs 
in analyzing three jet configurations that are close to being two jet configurations. 
On the other hand, the numerical method offers evident advantages in flexibility to 
modify the integrand. Since this method is quite new, one cannot yet say for what 
problems it might do better than the now standard ERT method. 

The numerical integration method exists as computer code with accompanying 
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technical notes and many of the basic ideas behind it have been described in two 
papers 0,^- In this talk, I briefly review the basics of the numerical integration 
method. Then, I display some graphs that illustrate the cancellation of singularities 
that occurs inside the integrand in the numerical method. Finally, I discuss some 
avenues for future research. 



2 Review of the numerical method 

Let us begin with a precise statement of the problem. We consider an infrared 
safe three-jet-like observable in e^e~ — > hadrons, such as a particular moment of the 
thrust distribution. The observable can be expanded in powers of a^/vr, 

a = ^aM, oc K/tt)" . (1) 

n 

The order contribution has the form 

1 f da^^^ 
= ^JdP^dP^dfM^'^'^'''^'^ 

+ ydp,dp,dp, Ss{puP,,ps) (2) 

+ — / dpidp2dp3dp4 ,^ ,^ ,^ 54(pi,P2,P3,P4). 

4! J apiap2dpsap4 

Here the dajf' are the order contributions to the parton level cross section, cal- 
culated with zero quark masses. Each contains momentum and energy conserving 
delta functions. The da^^^ include ultraviolet renormalization in the MS scheme. The 
functions S describe the measurable quantity to be calculated. We wish to calculate 
a "three-jet-like" quantity. That is, S2 = 0. The normalization is such that iS„ = 1 
for n = 2, 3, 4 would give the order perturbative contribution the the total cross 
section. There are, of course, infrared divergences associated with Eq. (^. For now, 
we may simply suppose that an infrared cutoff has been supplied. 

The measurement, as specified by the functions Sn, is to be infrared safe, as 
described in Ref. 0: the Sn are smooth, symmetric functions of the parton momenta 
and 

Sn+l{pi, . . . , Xpn, (1 - A)p„) = Sn{pi, ■ ■ ■ , Pn) (3) 

for < A < 1. That is, coUinear splittings and soft particles do not affect the 
measurement. 

It is convenient to calculate a quantity that is dimensionless. Let the functions 
Sn be dimensionless and eliminate the remaining dimensionality in the problem by 
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Figure 1: Two cuts of one of the Feynman diagrams that contribute to e^e hadrons. 



dividing by ctq, the total e^e cross section at the Born level. Let us also remove the 
factor of (as/Tr)^. Thus, we calculate 

X= , , (4) 

Let us now see how to set up the calculation of X in a convenient form. We note 
that X is a function of the cm. energy ^/s and the MS renormalization scale /x. We 
will choose /i to be proportional to ^/s■. fi = Af/yy^. Then X depends on Auv- But, 
because it is dimensionless, it is independent of ^/s. This allows us to write 

X= / dv^/l(v^) X(A^y,v^), (5) 

Jo 

where h is any function with 

d^s h{^s) = 1. (6) 







The quantity X can be expressed in terms of cut Feynman diagrams, as in Fig. |T[ 
The dots where the parton lines cross the cut represent the function iS„(pi, . . . ,Pn)- 
Each diagram is a three loop diagram, so we have integrations over loop momenta 
li, I2 and ^3. We first perform the energy integrations. For the graphs in which 
four parton lines cross the cut, there are four mass-shell delta functions These 
delta functions eliminate the three energy integrals over /°, 1 21 and /g as well as the 
integral over ^/s. For the graphs in which three parton lines cross the cut, we 
can eliminate the integration over ^/s and two of the integrals. One integral over 
the energy E in the virtual loop remains. We perform this integration by closing the 
integration contour in the lower half E plane. This gives a sum of terms obtained 
from the original integrand by some simple algebraic substitutions. Having performed 
the energy integrations, we are left with an integral of the form 

X = ^ / dh dh dh E 9{G, C- fi, h). (7) 

G '' C 
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Here there is a sum over graphs G (of which one is shown in Fig. |l|) and there is a 
sum over the possible cuts C of a given graph. The problem of calculating X is now 
set up in a convenient form for calculation. 

If we were using the Ellis- Ross- Terr ano method, we would put the sum over cuts 
outside of the integrals in Eq. (|^). For those cuts C that have three partons in the 
final state, there is a virtual loop. We can arrange that one of the loop momenta, say 
/i, goes around this virtual loop. The essence of the ERT method is to perform the 
integration over the virtual loop momentum analytically ahead of time. The integra- 
tion is often ultraviolet divergent, but the ultraviolet divergence is easily removed by 
a renormalization subtraction. The integration is also typically infrared divergent. 
This divergence is regulated by working in 3 — 2e space dimensions and then taking 
e ^ while dropping the 1/e" contributions (after proving that they cancel against 
other contributions). After the li integration has been performed analytically, the 
integrations over I2 and can be performed numerically. For the cuts C that have 
four partons in the final state, there are also infrared divergences. One uses either 
a 'phase space slicing' or a 'subtraction' procedure to get rid of these divergences, 
cancelling the pieces against the l/e" pieces from the virtual graphs. In the end, 
we are left with an integral / dli dl2 dl^ in exactly three space dimensions that can be 
performed numerically. 

In the numerical method, we keep the sum over cuts C inside the integrations. We 
take care of the ultraviolet divergences by simple renormalization subtractions on the 
integrand. We make certain deformations on the integration contours so as to keep 
away from poles of the form l/[Ep — Ej + ie]. Then the integrals are all convergent 
and we calculate them by Monte Carlo numerical integration. 

Let us now look at the contour deformation in a little more detail. We denote 
the momenta {h,l2,h} collectively by / whenever we do not need a more detailed 
description. Thus 



For cuts C that leave a virtual loop integration, there are singularities in the integrand 
of the form Ep — Ej + ie (or Ep — Ej — ie if the loop is in the complex conjugate 
amplitude to the right of the cut). Here Ep is the energy of the final state defined by 
the cut C and Ej is the energy of a possible intermediate state. These singularities do 
not create divergences. The Feynman rules provide us with the ie prescriptions that 
tell us what to do about the singularities: we should deform the integration contour 
into the complex / space so as to keep away from them. Thus we write our integral 
in the form 



Here in is a purely imaginary nine- dimensional vector that we add to the real nine- 
dimensional vector / to make a complex nine- dimensional vector. The imaginary part 





X = Y, dlY, J(G, C- 1) g(G, C- 1 + ik(G, C; /)). 



(9) 



G '' C 
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K depends on the real part /, so that when we integrate over /, the complex vector l+in 
lies on a surface, the integration contour, that is moved away from the real subspace. 
When we thus deform the contour, we supply a jacobian J = det{d{l + (See 
Ref. |] for details.) 

The amount of deformation k, depends on the graph G and, more significantly, the 
cut C. For cuts C that leave no virtual loop, each of the momenta li, I2, and I3 flows 
through the final state. For practical reasons, we want the final state momenta to be 
real. Thus we set k = for cuts C that leave no virtual loop. On the other hand, 
when the cut C does leave a virtual loop, we choose a non-zero k. We must, however, 
be careful. When k = there are singularities in g on certain surfaces that correspond 
to collinear parton momenta. These singularities cancel between g for one cut C and 
g for another. This cancellation would be destroyed if, for / approaching the collinear 
singularity, k, = for one of these cuts but not for the other. For this reason, we 
insist that for all cuts C, k — as / approaches one of the collinear singularities. The 
details can be found in Ref. [Q. 

Much has been left out in this brief overview, but we should now have enough 
background to see how the method works. 



3 Example 

I present here a simple example, taken from Ref. ||^. Instead of working with 
QCD at three loops with many graphs, let's work with one graph for (p^ theory at 
two loops, as shown in Fig. 0. 




Figure 2: Sample graph in (p"^ theory. 



This graph has four final state cuts, as shown in Fig. |^. We will fix the incoming 
momentum q and integrate over the incoming energy For a measurement function, 
we take S{p) = J2 \PT,i\, where px is the part of p orthogonal to q. We make a choice 
of contour deformations and of the density p of Monte Carlo integration points as 
described in |^ . Then we can plot the integrand / divided by the density of points p 
versus the loop momentum. In a Monte Carlo integration, large f/p corresponds to 
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Figure 3: Cuts of the sample graph. 



large fluctuations, so //p should never be too large. In the two figures that follow, 
I plot f /p versus the momentum in the left hand loop. Specifically, using the kn 
defined in Fig. H, I plot//p versus / = /c2 at fixed /C4 for / in the {/C4, q\ plane. 




Figure 4: Integrand divided by the density of points for the three parton cuts. The collinear 
singularities are visible. 

In Fig. ^, I show f / p summed over the two cut Feynman graphs that have three 
partons in the final state, leaving no virtual loop. Evidently, there are singularities. 
There is a soft parton singularity (at / = 0) that I have cut out of the diagram 
and there are collinear parton singularities that are visible in the picture. In the 
Ellis-Ross- Terrano method, these cut graphs would be calculated using a numerical 
integration. But first a cutoff or other method for eliminating the singularities would 
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be needed to eliminate the singular region. The two cuts that leave virtual subgraphs 
also lead to singularities along the coUinear lines in the space of the loop momentum. 
I omit displaying a graph of //p for these two cut Feynman graphs because the 
result simply looks like an upside down version of Fig. ^. In the Ellis-Ross- Terrano 
method, one takes care of the singularities in the virtual loop by integrating in 3 — 2e 
space dimensions. In the numerical method, one combines the integrands for all of the 
cuts. Then the coUinear singularities disappear, while the soft singularity is weakened 
enough that it can be eliminated in //p by building a suitable singularity into p. As 
suggested by the title for this talk, the cancellation of singularities between real and 
virtual graphs happens by itself because it is built into the Feynman rules. The result 
for //p summed over all four cuts is shown in Fig. ^. The coUinear singularities are 
gone, while the soft parton singularity in / has been weakened enough that it is 
cancelled by a corresponding singularity in p. Thus a Monte Carlo integration of / 
using a density of integration points p can converge nicely because f /pis not singular. 

What remains visible in Fig. ^ is a ridge in //p for / lying on the ellipsoidal surface 
defined by |fci| + j/cal = |fc4| + l^sl, where the intermediate state energy in the virtual 
graphs matches the final state energy. This ridge is related to an energy denominator 
factor l/lEp — Ei+ie] in old fashioned perturbation theory. The numerical integration 
method has taken advantage of the ie prescription in the Feynman rules to deform 
the integration contour and avoid the singularity. 




Figure 5: Integrand divided by the density of points for all cuts together. The coUinear 
singularities disappear while the soft parton singularity in / is weakened so that it can be 
cancelled by a singularity in p. 
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4 Prospects 



There are a number of promising areas for further research along these hues. 

The current program, beowulf [|], does e"*" + e~ — >• 3 jets at next-to-leading or- 
der. The partons are all massless. With some modifications, the partons could have 
masses. Then one could include massive quarks and one could extend the theory to 
the complete Standard Model with its massive vector bosons. Furthermore, one could 
add supersymmetry interactions. 

The current program is confined to processes with no hadrons in the initial state. 
Presumably the same idea can be applied to processes with initial state hadrons, 
that is electron-proton collisions and proton-proton or proton-antiproton collisions. 
Again, one should be able to make the particles massive so that one can extend the 
calculations to the complete Standard Model and supersymmetry. 

It should also be possible to have more final state partons. That is, one could 
attempt to calculate e"*" -|- e~ — > 4 jets or p -|- p ^ 3 jets at next-to-leading order. 

The challenge of the legendary hero Beowulf was to kill the monster Grendel. The 
monsters listed above are already dead or at least gravely injured. In particular, all 
that beowulf can do could already have been accomplished by the program of Kunszt 
and Nason eleven years ago [Q. However, the challenge of calculating e"'" + e~ — 3 jets 
at next-to- next-to-leading order remains unmet, and it may be that a completely 
numerical attack would be successful. 

A less difficult goal is to use the flexibility inherent in the numerical method to go 
beyond fixed order perturbation theory. For instance, one could use running couplings 
inside the next-to-leading order graphs as a method for investigating power suppressed 
( "renormalon" ) contributions to the theory. More importantly, one could put a next- 
to-leading order calculation inside a parton shower event generator (or attach parton 
showers to the outside of the next-to-leading order calculation) in order to have a full 
parton shower event generator that is correct at next-to-leading order for three jet 
quantities in e~^e~ annihilation. This is, of course, not quite trivial 0. As far as I can 
see, the first step is to convert the current algorithms so that they operate in Coulomb 
gauge instead of Feynman gauge. In this way, the partons propagating into the final 
state have physical polarizations only. Then these physically polarized partons can 
split many times to make parton showers. One simply has to avoid counting the same 
splittings twice. 
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